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Abstract 

Within any given tissue, gene expression levels can vary extensively among indi¬ 
viduals. Such heterogeneity can be caused by genetic and epigenetic variability and 
may contribute to disease. The abundance of experimental data now enables the 
identification of features of gene expression profiles that are shared across tissues and 
those that are tissue-specific. While most current research is concerned with charac¬ 
terizing differential expression by comparing mean expression profiles across tissues, 
it is believed that a significant difference in a gene expression’s variance across tis¬ 
sues may also be associated with molecular mechanisms that are important for tissue 
development and function. 

We propose a sparse multi-view matrix factorization (sMVMF) algorithm to jointly 
analyse gene expression measurements in multiple tissues, where each tissue provides 
a different ‘view’ of the underlying organism. The proposed methodology can be 
interpreted as an extension of principal component analysis in that it provides the 
means to decompose the total sample variance in each tissue into the sum of two 
components: one capturing the variance that is shared across tissues and one iso¬ 
lating the tissue-specific variances. sMVMF has been used to jointly model mRNA 
expression profiles in three tissues obtained from a large and well-phenotyped twins 
cohort, TwinsUK. Using sMVMF, we are able to prioritize genes based on whether 
their variation patterns are specific to each tissue. Furthermore, using DNA methy- 
lation profiles available, we provide supporting evidence that adipose-specific gene 
expression patterns may be driven by epigenetic effects. 
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1 Introduction 


RNA abundance, as the results of active gene expression, affects cell differentiation 
and tissue development (Coulon et al, 2013). As such, it provides a snapshot of the 
undergoing biological process within certain cells or a tissue. Except for house-keeping 
genes, the expressions of a large number of genes vary from tissue to tissue, and some 
may only be expressed in a particular tissue or a certain cell type (Xia et al., 2007). 
The regulation of tissue-specific expression is a complex process in which a gene’s 
enhancer plays a key role regulating gene expressions via DNA methylation (Ong and 
Corces, 2011). Genes displaying tissue-specific expressions are widely associated with 
cell type diversity and tissue development (Reik, 2007), and aberrant tissue-specific 
expressions have been associated with diseases that originated in the underlying tissue 
(van’t Veer et al., 2002; Lage et al., 2008). Distinguishing tissue-specific expressions 
from expression patterns prevalent in all tissues holds the promise to enhance fun¬ 
damental understanding of the universality and specialization of molecular biological 
mechanisms, and potentially suggest candidate genes that may regulate traits of in¬ 
terest (Xia et al, 2007). As collecting genome-wide transcriptomic profiles from many 
different tissues of a given individual is becoming more affordable, large population- 
based studies are being carried out to compare gene expression patterns across human 
tissues (Liu et al, 2008; Yang et al., 2011). 

A common approach to detecting tissue-specific expressions consists of comparing 
the mean expression levels of individual genes across tissues. This can be accom¬ 
plished using standard univariate test statistics. For instance, Wu et al. (2014) used 
the two-sample Z-test to compare non-coding RNA expressions in three embryonic 
mouse tissues: they reported approximately 80% of validated in vivo enhancers ex¬ 
hibited tissue-specific RNA expression that correlated with tissue-specific enhancer 
activity. Yang et al. (2011) applied a modified version of Tukey’s range test (Tukey, 
1949), a test statistic based on the standardised mean difference between two groups, 
to compare expression levels of 127 human tissues, and results of this study are pub¬ 
licly available in the VeryGene database. A related database, TiGER (Liu et al., 
2008), has also been created by comparing expression sequence tags (EST) in 30 
human tissues using a binomial test on EST counts. Both VeryGene and TiGER con¬ 
tain up-to-date annotated lists of tissue-specific gene expressions, which generated 
hypotheses for studies in the area of pathogenic mechanism, diagnosis, and therapeu¬ 
tic research (Wu et al, 2009). 

More recent studies have gone beyond the single-gene comparison and aimed at 
extracting multivariate patterns of differential gene expression across tissues. Xiao 
et al. (2014) applied the higher-order generalised singular value decomposition (HO- 
GSVD) method proposed by Ponnapalli et al. (2011) and compared co-expression 
networks from multiple tissues. This technique is able to highlight co-expression pat¬ 
terns that are equally significant in all tissues or exclusively significant in a particular 
tissue. The rationale for a multivariate approach is that when a gene regulator is 
switched on, it can raise the expression level of all its downstream genes in specific 
tissues. Hence a multi-gene analysis may be a more powerful approach. 

While most studies explore the differences in the mean of expression, the sam- 
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pie variance is another interesting feature to consider. Traditionally, comparison of 
expression variances has been carried out in case-control studies (Mar et al., 2011). 
Using an F-test, significantly high or low gene expression variance has been observed 
in many disease populations including lung adenocarcinoma and colerectal cancer, 
whereas the difference in mean expression levels was not found significant between 
cases and controls (Ho et al., 2008). In a tissue-related study, Cheung et al. (2003) 
carried out a genome-wide assessment of gene expressions in human lymphoblastoid 
cells. Using an F-test, the authors showed that high-variance genes were mostly 
associated with functions such as cytoskeleton, protein modification and transport, 
whereas low-variance genes were mostly associated with signal transduction and cell 
death/proliferation. 

In this work we introduce a novel multivariate methodology that can detect pat¬ 
terns of differential variance across tissues. We regard the gene expression profiles in 
each tissue as providing a different “view” of the underlying organism and propose 
an approach to carry out such a multi-view analysis. Our objective is to identify 
genes that jointly explain the same amount of sample variance in all tissues - the 
"shared" variance - and genes that explain substantially higher variances in each 
specific tissue separately - the "tissue-specific" variances - while the shared variance 
has been accounted for. During this process we impose a constraint that the factors 
driving shared and tissue-specific variability must be uncorrelated so that the total 
sample variance can be decomposed into the two corresponding components. The 
proposed methodology, called sparse multi-view matrix factorisation (sMVMF), can 
be interpreted as an extension of principal component analysis (PCA), which is tra¬ 
ditionally used to identify a handful of latent factors explaining a large portion of 
sample variance separately in each tissue. 

The rest of this paper is organised as follows. The sMVMF methodology is pre¬ 
sented in Section 2, where we also discuss connections with a traditional PCA and 
derive the parameter estimation algorithm. In Sectionv 3 we demonstrate the main 
feature of the proposed method on simulated data, and report on comparison with al¬ 
ternative univariate and multivariate approaches. In Section 4 we apply the sMVMF 
to compare mRNA expressions in three tissues obtained from a large twin population, 
the TwinsUK cohort. We conclude in Section 5 with a discussion. 


2 Methods 

2.1 Sparse multi-view matrix factorisation 

We assume to have collected p gene expression measurements for M different tissues. 
Ideally the data for all tissues should be derived from the same underlying random 
sample (as in our application, Section 4) in order to remove sources of biological 
variability that can potentially induce differences in gene expression profiles across 
tissues. In practice, however, cross-tissue experiments rarely collect samples from the 
same set of subjects or may fail quality control. In our setting therefore we assume M 
different random samples, each one contributing a different tissue dataset. The m th 
dataset consists of n m subjects, and the expression profiles are arranged in an n m x p 
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matrix. All matrices are collected in X = {AlW, X^ 2 \ ..., where the super¬ 

scripts refer to tissue indices. For each X^ m \ we subtract the column mean from each 
column such that each diagonal entry of the scaled gram matrix, -X(x( m ^) T X( m \ is 
proportional to the sample variance of the corresponding variable, and the trace is the 
total sample variance. We aim to identify genes that jointly explain a large amount of 
sample expression variances in all tissues and genes that explain substantially higher 
variances in a specific tissue. Our strategy involves approximating each ^== X by 
the sum of a shared variance component and a tissue-specific component: 
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m r 


~ g( m ) _|_ (i) 

shared variance component tissue-specific variance component 


for m = 1,2,..., M, where 1 / \Jn m is a scaling factor such that the trace of the gram 
matrix of the left-hand-side equals the sample variance. These components are defined 
so as to yield the following properties: 

(a) The rank of S and T are both much smaller than nrin(n m ,p) so that the 
two components provide insights into the intrinsic structure of the data while 
discarding redundant information. 


(b) 


(c) 


The variation patterns captured by shared component are uncorrelated to the 
variation patterns captured by tissue-specific component. As a consequence of 
this, the total variance explained by S and T^ altogether equals the sum 
of the variance explained by each individual component. 


The shared component explains the same amount of variance of each gene ex¬ 
pression in all tissues. As such, the difference in expression variance between 
tissues is exclusively captured in tissue-specific variance component. 

We start by proposing a factorisation of both S and T^ which, by imposing 
certain constraints, will satisfy the above properties. Suppose rank(S'^ m ^) = d and 

rank(T(”b) = rf where d, r « min (n m ,p) following property (a). For a given r, 
T (m) 

can be expressed as the product of an n m x r full rank matrix W^ and the 
transpose of a p X r full rank matrix V^ m \ that is: 


= jyWjyHjT = y( 


3 = 1 


E 7 ] 

3 = 1 


\j] 


( 2 ) 


where the superscript T denotes matrix transpose, and the subscript j denotes the 
j th column of the corresponding matrix. Each 

:= wj m \vj m ' > ) T 

has the same dimension as T*™-) and is composed of a tissue-specihc latent factor 
(LF). A LF is an unobservable variable assumed to control the patterns of observed 
variables and hence may provide insights into the intrinsic mechanism that drives the 
difference of expression variability between tissues. The matrix factorisation in (2) is 
not unique, since for any r X r non-singular square matrix R, T( m ' ) = W(yi m )^ T = 
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(W^R)^^ 1 (V ( ' m ' l ) T ) = (y( m ')) T . We introduce an orthogonal constraint 

( W^ m ^) T W ^ m ) = I r so that the matrix factorisation is unique subject to an isometric 
transformation. Similarly, we can factorise the shared component as: 


S {m) = jj(m )(y*)T = {V£) T = ^ S \ 


: (m) 

[k] 


(3) 


k= 1 


fc=l 


where U^ is orthogonal and V* is tissue-independent which we shall explain. Each 
Srjq has the same dimension as S^ and is composed of one shared variability LF. 
The resulting multi-view matrix factorisation (MVMF) then is: 

—-—xM ss U^ m \V*) T + W^ m) (V (m) ) T (4) 

The matrix factorisations (2) and (3) are intimately related to the singular value 
decomposition (SVD) of S and T^ m \ Specifically, U and W are analogous 
to the matrix of left singular vectors and also the principal components (PCs) in a 
standard PCA. They represent gene expression patterns in a low-dimensional space 
where each dimension is derived from the original gene expression measurements such 
that the maximal amount of variance is explained. We shall refer the columns of U^ 
and W as the principal projections (PPJ). ( V*) T and are analogous to 

the product of the diagonal matrix of eigenvalues and the matrix of right singular 
vectors. Since the singular values determine the amount of variance explained and 
the right singular vectors correspond to the loadings in the PCA which quantifies 
the importance of the genes to the expression variance explained, using the same 
matrix V* for all tissues in the shared component results in the same amount of 
shared variability explained for each gene expression probe, such that property (c) is 
satisfied. We shall refer to matrices V* and f/l m ) as transformation matrices. 

A sufficient condition to satisfy property (b) is: 

= 0 dxr (5) 

This constraint, in addition to the orthogonality of U^ and W^ m \ results in the ( d+ 
r ) PPJs represented by being pairwise orthogonal, which is analogous 

to the standard PCA where the PCs are orthogonal. Intuitively, this means for each 
tissue the LFs driving shared and tissue-specific variability are uncorrelated. The 
amount of variance explained in tissue m, can be computed as (subject to a 
constant factor): 

d* m = Tr {(S (m) ) T S (m) + ( T M) r T+ 2{S im) ) T T {m] } (6) 

where Tr denotes the matrix trace. Recalling that S4 m ) = and ([/( m )) T [/"( m ) 

Id, the amount of shared variance explained is: 

a * = Tr{(S( m )) T S( m )} = Ty{V*{V*) t } (7) 


Likewise, recalling that p( m ) = W <ym \V^ m ^) T and = I r , the amount 

of tissue-specific variance explained is: 

a m = Tr {(T (m )) T T( m )} = Tr{V^ m) (R (m) ) T } (8) 
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Making the same substitutions into (6), we obtain: 


a s m = Tv{V*(V*) t + V {m \V^ m) ) T + 2 V*(U {m) ) T W {m \V^ m) ) T } 

Substituting (5) into the above equation, we reach: 

a s m = Tr{V*(V*) T + V^ m \V^) T } = + °m 


(9) 


which satisfies (b). 

2.2 Sparsity constraints and estimation 

The factorisation (4) is obtained by minimising the squared error. This amounts to 
minimising the loss function: 

M 

£ = II —— x(m) - U {m) (V*) T - W {m) [V {m) ) T ||3r (10) 

m=l v 

where ||.||jr refers to the Frobenius norm, subject to the following orthogonality con¬ 
straints: 


( [ /M)T [/ (m) = j ( W (m))T w (m) = j = Q (H) 

For fixed U^ m \V*) T , the optimal T= W^(V^ m ^) T is a low-rank approximation 
of II — S'^ m ^||?r, where each rank sequentially captures the maximal variance 

remained in each data matrix after removing the shared variability. Likewise, for fixed 
, each rank of the optimal S( m ' ) = U^ m \V*) T sequentially captures 
the maximal variance remained across all tissues after removing the tissue-specific 
variance. 

In transcriptomics studies, it is widely believed that the differences in gene expres¬ 
sions between cell and tissue types are largely determined by transcripts derived from 
a small number of tissue-specific genes (Jongeneel et al., 2005). Therefore it seems 
reasonable that in our application of multi-tissue comparison of gene expressions, for 
each PPJ, the corresponding column in the transformation matrix should feature a 
limited number of non-zero entries. In such a scenario, a sparse representation will 
not only generate more reliable statistical models by excluding noise features, but also 
offer more biological insight into the underlying cellular mechanism (Ma and Huang, 
2008). 

In the context of MVMF, we induce sparse estimates of V* and y( m ) by adding 
penalty terms to the loss function l (U, W, V*, V) as in (10). Specifically, we minimise: 

M 

£ {U,W,V*,V) + 2 ■ M ■ ||H*Ali + 2^ ||H (m) A (m) ||i (12) 

771—1 

where || ||i denotes the £\ norm. A* and A( m ) are d X d and r x r diagonal matrices, 
respectively. In both matrices, the k th diagonal entry is a non-negative regularisation 
parameter for the k th column of the corresponding transformation matrix, and the 
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k th column tends to have more zero entries as the k th diagonal entry increases. In 
practice, a parsimonious parametrisation may be employed where A* = \\Id and 
AW = A 2 J r for m = 1 so that the number of parameters to be specified is 

greatly reduced. Alternatively, A* and A( m ) may be set such that a specified number 
of variables are selected in each column of V* and V^ m \ 

The optimisation problem (12) with constraints (11) is not jointly convex in U^ m \ 

w {m)^ y(m) ) and V * f or 

m = 1,2(for instance the orthogonality constraints 
are non-convex in nature), hence gradient descent algorithms will suffer from multiple 
local minima (Gorski et al ., 2007). We propose to solve the optimisation problem by 
alternately minimising with respect to one parameter in U^ m \W^ m \V* , while 

fixing all remaining parameters, and repeating this procedure until the algorithm con¬ 
verges numerically. The minimisation problem with respect to V* or Vd m ) alone is 
strictly convex, hence in these steps a coordinate descent algorithm (CDA) is guar¬ 
anteed to converge to the global minimum (Friedman et al. , 2007). CDA iteratively 
update the parameter vector by cyclically updating one component of the vector at a 
time, until convergence. On the other hand, the minimisation problem with respect to 
W or [/( m ) is not convex. For fixed V* and V^ m \ the estimates of W( m ' ) and U 
that minimise (12) can be jointly computed via a closed form solution. Assuming we 
have obtained initial estimates of V* and we cyclically update the parameters 

in the following order: 

(t/( m) , W (m) ) ->• ->• 

Here U {m "> and W^ are jointly estimated in the first step, and in the subsequent 
steps V ( m ^ and V* are updated separately, while keeping the previous estimates fixed. 
A detailed explanation of how each update is performed is in order. 

First we reformulate the estimation problem as follows: we bind the columns of 
U and W and define the n m x ( d+r ) augmented matrix: f7^ m) = [U^ , 
we then bind the columns of V* and I/( m ) and define the p x (d + r) matrix: t/! m ) = 
[V* , V (m) ]. As such: 

M 

£ (U,W,V*,V {m) ) = V II-- C/ (m) (H (m) ) T || 2 F 


and the constraints in (11) can be combined into: 

(C/Mft/M = I d+r 

Fixing V^ m \ the estimate of U can be obtained by the reduced-rank Procrustes 
rotation procedure which seeks the optimum rotation of X such that the error 
||^=AT( m ) — F r ( m )(V r ( m ))' i ||2_ j s minimal. For a proof of this, see (Zou et al, 2006). 

We obtain the SVD of as PQR T , and compute the estimate of U^ 

by: [/M = PR t . 

Next, we fix U^ m \ W^ m \ and V* while minimising (12) with respect to . 
For each fixed m, varying I/( m ) only changes the objective function via the summand 
indexed (m). Hence it is sufficient to minimise: 


1 

_ z _ x - U^ m \V*) T 

\Jr+n 


Wi™) (y( m )) T || 2 _ + 2||H( m )A (m) ||i 


(13) 
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This function is strictly convex in and the CDA is guaranteed to converge to 

the global minimum. We drop the superscript (m) in the following derivation for 
convenience and denote the j th column of the matrix V by Vj. In each iteration, the 
estimate of Vj is found by equating the first derivative of (13) with respect to Vj to 
zero. Hence: 


- 2 ( 



X - UV* 


WV T ) T Wj + 2Aj ■ V(\Vj\) = 0, 


where V is the gradient operator. Substitute (11) and rearrange to give: 

1 


Vj = 


/1l r 


-X Wj — Aj ■ V(|V}|) 


We dehne the sign function <r(y) which equals 1 if y > 0, —1 if y < 0, and 0 if y = 0. 
First note the derivative of the function |y| is cr(y) if y ^ 0 and a real number in the 
interval ( — 1,1) otherwise. Rearrange the previous equation to obtain the updated 
estimate in each iteration: 

v} m) = S A (m) ((^=X (m ^ T wj m) ) (14) 

3 y Tim 


where S\(y) is a soft-thresholding function on vector y with non-negative parameter 
A such that S\(y) = <r(y) ■ max{|y| — A, 0}, and is the j th diagonal entry of A^ m ^. 

In the third step, we fix the estimates of U^ m \ W^ m \ and V^ and minimise (12) 
with respect to V* . The objective function becomes: 


1 + 2 • M ■ || V*A* ||i (15) 

where ^ is defined in (10). As in the second step, we use a CDA in each iteration and 
the updated estimate of V* is found by equating the first derivative of (15) to zero. 
Specifically: 


- 2 Em=i - U^V* - W ( ' m \V ( ' m ' > ) T ] T Ui^j 

+ 2 • M • A* • V(|TA*|) = 0, 

where A* is the i th diagonal entry of A*. Applying (11), this can be re-arranged into: 


M ■ V* 


M 

\/ 1 l m, 


Tjj{rn) 


m=1 


M-A*-V(\V*\), 


Using the soft-thresholding and the sign functions, the updated estimate in each 
iteration can be re-written as: 


V* 

v l 



(16) 


The cyclic CDA requires initial estimates of V* and V^ m \ which are obtained as 
follows. First we set an initial value to V*, which explains as much variance in all 








datasets in X as possible. This amounts to a PCA on the (X)m=i n m) x P ma tri x V 
obtained by binding the rows of m = 1, M. We compute the truncated 

V 

SVD of X and obtain X = UDB T where D contains the d largest eigenvalues of 
X T X. The initial estimate of V* is then defined as: 

0n T = jjDB T , (17) 

and is defined by the corresponding rows of U in the SVD. For the tissue-specific 
transformation matrices W m ), we compute the SVD of the residuals after removing 
the shared variance component from }— X^ m \ which gives: i— X^ — = 

Y Tim v n m 

W( m )i?( m )(Q( m )) T . The initial estimate of is defined as: 

= i?W(QW) T . (18) 

A summary of the estimation procedure is given in Algorithm 1. 


Algorithm 1 sMVMF estimation algorithm 

Input: data X; parameters d, r, A*” 1 ), A* for m = 1, 2,..., M. 

Output: U^ m \ for m = 1,2, and V*. 

1: Get initial estimates of V for m = 1, 2,..., M, and V* as in (18) and (17). 
2: while not convergent do: 

3: Apply SVD: ^X (m )V (m ) = PQR T , and set U (m) = Pi? T . 

4: Use CDA to estimate V^ according to (14). 

5: Use CDA to estimate V* according to corollary (16). 


2.3 Parameter selection 

The sMVMF contains two sets of parameters: the tissue-specific sparsity parameters 
A*, and the (d,r) pair. Both d and r balance model complexity and the amount 
of variance explained. We select the smallest possible values of d and r such that a 
prescribed proportion of variance is explained. For a fixed (d, r) pair, the sparsity 
parameters can be optimised using a cross-validation procedure, which identifies the 
best combination from a grid of candidate values so that the amount of variance 
explained is maximised on the testing data for the chosen ( d,r). However, in high¬ 
dimensional settings, cross-validation procedures such as this one tend to favour over¬ 
complex models which may include noise variables (Biihlmann and van de Geer, 2011). 
Instead we propose using the “stability selection” procedure which is particularly 
effective in improving variable selection accuracy and reducing the number of false 
positives in high-dimensional settings (Meinshausen and Biihlmann, 2010). Given 
parameters d = do and r = ro in sMVMF, variables can be ranked according to their 
importance in explaining shared and tissue-specific variances by applying a stability 
selection procedure as follows: 
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1. Randomly extract half of the n m samples from each X without replacement 
and denote the resulting data matrix xj m \ for m = 1, M. In the case where 
each A"h n ) consists of the same subjects, it may be preferable to draw the same 
samples from all datasets. 

2. Fit sMVMF on Xg m \ rri = 1, M, where A* and A( m ) are chosen such that a 
prescribed number of variables are selected from each column of V* and V^ m \ 
m = 1, M. 

3. Record the variables that are selected in V* up to and including d = do and in 
y( m ) U p an d including r = ro, m = 1, 

4. Repeat steps 1 to 3 IV times, where N is at least 1000. 

5. Compute the empirical selection probabilities for each variable in V* and V^ m \ 

m = 1 Then rank the variables in each list according to the selection 

probabilities. 

Note in step 2, the number of variables to be selected in each column of V* and 
y(»), m = 1 is a regularisation parameter. Nevertheless, Meinshausen and 

Biihlmann (2010) showed the variable rankings, especially the top ranking variables, 
were insensitive to the choice of these parameters which regularised the level of spar¬ 
sity. In practice, the number of variables selected in each column of V* and V^ m \ 
rn = 1 , ...,M, is randomly picked and kept small. In the TwinsUK study, it is chosen 
to be 100 since we would only be interested in the top few hundred probes which 
drove the shared and tissue-specific variability respectively. 


3 Illustration with simulated data 

In this section we present simulation studies to characterise how the sMVMF method 
is able to distinguish between shared and tissue-specific variance. We simulate shared 
and tissue-specific variance patterns as illustrated by the middle and right panels in 
Figure 1. We then test whether sMVMF correctly decomposes the total sample vari¬ 
ance (left panel) whilst detecting variables contributing to the non-random variabil¬ 
ity within each variance component. We also compare sMVMF with two alternative 
methods: standard PCA and Levene’s test (Gastwirth et al . , 2009) of the equality of 
variance between population groups. 

3.1 Simulation setting 

Our simulation study consists of 1000 independent experiments. In each experiment 
we simulate 3 data matrices or datasets (tissues) of dimension n = 100 (samples) and 
p = 500 (genes). Each simulated data matrix X is obtained via: 

X M = y( m ) _|_ ^( m ) _|_ x( m ) ^ 

where yl m ) j s a component designed to control the shared variance, Z^ is introduced 
to control the tissue-specific variance, and E is a random error. They are all 
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n x p random matrices. Since we ultimately wish to test whether our method is able 
to distinguish between signal and noise variables, we assume that only the first 30 
variables carry the signal, whereas the remaining 470 only introduce noise. 

We suppose that the shared variability is controlled by the activation of 3 latent 
factors, each regulating the variance of a different block of variables. To this end, 
we further group the 30 signal variables into three blocks of 10 normally distributed 
random variables each (numbered 1 —10,11 — 20, and 21 — 30), as illustrated in Figure 
2 (A). We design the simulations so that each of the first 30 variables in Y has the 
same variance in different datasets; moreover, the variance decreases while moving 
from the first to the third block. Further details and simulation parameters are given 
in Appendix, Section A. This procedure generates shared variance patterns that look 
like those reported in the middle panel of Figure 1. 

The variables in Z are also assumed to be normally distributed. They are gen¬ 
erated such that exactly 10 of them have the largest variance across datasets. The 
resulting "mosaic" structure of the simulated variance patterns is illustrated in right 
panel of Figure 1. The data matrices yb ra ) and ^( m ) are generated such that the total 
non-random sample variance of each variable in a tissue equals the sum of its shared 
and tissue-specific variances, which is also illustrated in Figure 1. The random error 
term is generated from independent and identical normal distributions with zero 
mean and noise of for all variables in all datasets. We perform simulations on two 
settings: in setting I of = 1 and in setting II (j'f = 4. As a result of this simulation 
design, we are able to characterise the true underlying architecture that explains the 
total sample variance. 

3.2 Simulation results 

The data generated in each experiment was analysed by fitting the sMVMF algorithm. 
To focus on the ability of the model to disentangle the true sources of variability, we 
take d = 3 and r = 1, which equal the true number of shared and tissue-specific LFs 
used to generate the data. The regularisation parameters A* and A( m ) are tuned such 
that each PPJ consists of 10 variables, the true number of signal variables. 

For comparison, we propose two additional approaches that are able to identify 
variables featuring dataset-specific sample variances, although they do not attempt 
to model the shared variance. The first method consists of carrying out a separate 
PCA on each dataset; for each PCA/dataset, we then select the 10 variables having 
the largest loadings in the first principal component. The second method consists of 
applying a standard Levene’s test of equality of population variances independently 
for each variable, which is then followed by a Bonferroni adjustment to control the 
family-wise error rate; if a test rejects the null hypothesis at the 5% significance level, 
we select the variable having the largest sample variance amongst the three datasets. 

By averaging across 1000 experiments, we are able to estimate the probability 
that each one of the 30 signal variables is selected by each one of the three compet¬ 
ing methods. The heatmaps (A)-(C) in Figure 3 visually represent these selection 
probabilities for simulation setting I. Here sMVMF perfectly identifies the variables 
that introduce dataset-specific variability. The results obtained using Levene’s tests 
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Tissue 1 
Tissue 2 
Tissue 3 

Figure 1: Simulated patterns of sample variance: the total, non-random, sample vari¬ 
ance of 30 signal-carrying random variables is generated so that it can be decomposed 
into the sum of shared and tissue-specific components. Rows correspond to tissues 
(datasets) and columns correspond to 30 variables. Brighter colours represent large 
variance and darker colours represent low variance. Although by construction the un¬ 
derlying shared and tissue-specific variances have very different patterns, sMVMF is 
able to discriminate between them. 

True and estimated latent structure controlling shared variance 
(A) shared LF structure (B) sMVMF (C) Stacked-PCA 

Top to I ■ HHH 

Top 20 T Til 

Top 30 

Active Inactive o i 

Figure 2: Each latent factor (LF) is only active in a block of 10 signal- carry¬ 
ing variables, and controls the amount of variance of those variables that is shared 
amongst datasets. The (A) panel shows the true latent structure used to generate the 
data. Panels (B) and (C) show the estimated probabilities that each variable has been 
selected as signal-carrier using sMVMF and a stacked-PCA approach, respectively. 
sMVMF accurately captures the true shared LF structure whereas stacked-PCA tends 
to identify variables with large variance but fails to identify the LF structure. 



Decomposition of simulated total variance 
Total variance Shared variance Tissue-specific variance 



are somewhat similar, except for some variables in the first block (indexed 3 — 8) 
and second block (indexed 14 — 17). By reference to the middle panel of Figure 1, it 
can be noted that these variables are precisely those featuring large shared variabil¬ 
ity by construction. On the other hand, the PCA-based approach performs poorly 
because it can only select variables that contribute to explaining the total sample 
variance, but is unable to capture dataset-specific patterns. This example is meant 
to illustrate the limitations of both univariate and multivariate approaches that do 
not explicitly account for factors driving shared and dataset-specific effects. sMVMF 
has been designed to address exactly these limitations. 

Both Levene’s test and the individual-PCA approach are not designed to capture 
shared variance patterns. As a way of direct comparison with sMVMF we therefore 
propose an alternative PCA-based approach that has the potential to identify vari¬ 
ables associated to the direction of largest variance across all three datasets. This 
method consists of performing a single PCA on a “stacked” matrix of dimension 
{Ain) x p containing measurements collected from all three datasets, and obtained 
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Selected variables controlling tissue-specific variance 

Tissue 1 
Tissue 2 
Tissue 3 


Figure 3: Three different methods - sMVMF, Levene’s test and PCA - are used 
to detect random variables whose variance pattern is dataset-specific. Each heatmap 
represents the selection probabilities estimated by each method: (A) sMVMF produces 
patterns that closely match the true tissue-specific variances shown in the right panel 
of Figure 1; (B) Levene’s test performs well for variables those variance is mostly 
driven by tissue-specific factors, but fails to detect those variables having a strong 
shared-variance component; (C) The PCA-based method cannot distinguish between 
shared and tissue-specific variability, and fails to recover the true pattern. 


(A) sMVMF (B) Levene's test (C) PCA 



o,2 sj ce o« 


by coalescing the rows of the three individual data matrices. By varying the cutoff 
value for thresholding the loadings of the first PC, we are able to select the top 10, 
20, and 30 variables. We shall refer to this approach as stacked-PCA. 

Results produced by sMVMF and stacked-PCA are summarised by the heatmaps 
(B) and (C) in Figure 2, and can be directly compared to the true simulated patterns 
in (A). As expected, stacked-PCA tends to select variables having large total sample 
variances, whereas sMVMF can identify variables affected by each shared LF which 
jointly explain a large amount of variance. This example shows that sMVMF is able to 
identify the variables associated to the latent factors controlling the shared variance. 

We also carried out a simulation, based upon the same setting, with smaller signal- 
to-noise ratio, i.e. by sampling the random error terms in E from independent nor¬ 
mal distributions having larger variance. The results were very similar to the previous 
setting, except that Levene’s test was hardly able to identify any tissue-specific genes. 
The heatmaps summarising model performances are given in Appendix, Section B. 


4 Application to the TwinsUK cohort 

4.1 Data preparation 

TwinsUK is one of the most deeply phenotyped and well-characterised adult twin 
cohort in the world (Moayyeri et al., 2013). It has been widely used in studying the 
genetic basis of aging procession as well as complex diseases (Codd et al., 2013). More 
importantly, it contains a broad range of ‘onrics’ data including genomic, epigenomic 
and transcriptomic profiles amongst others (Bell et al, 2012). In this study, we 
focus on comparing the variance of mRNA expressions in adipose (subcutaneous 
fat), lymphoblastoid cell lines (LCL), and skin tissues. The microarray data used 
in this study were obtained from the Multiple Tissue Human Expression Resource 
(Nica et al., 2011), with participants being recruited from the TwinsUK registry. 
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Peripheral blood samples were artificially transformed from mature blood cells by 
infecting them with the Epstein-Barr virus (Glass et al., 2013). All tissue samples 
were collected from 856 female Caucasian twins (154 monozygotic twin pairs, 232 
dizygotic twin pairs and 84 singletons) aged between 39 and 85 years old (mean 
62 years). Genome-wide expression profiling was performed using Illumina Human 
HT-12 V3 BeadChips, which included 48,804 probes. Log2-transformed expression 
signals were normalized per tissue using quantile normalization of the replicates of 
each individual followed by quantile normalization across all individuals, as described 
in Nica et al. (2011). In addition, we also had access to 450K methylation data of the 
same adipose biopsies profiled using Inhnium HumanMethylation 450K BeadChip Kit 
(Wolber et al., 2014). We only retained probes whose expression levels were measured 
in all three tissues, and removed subjects comprising unmeasured expressions in any 
tissue. Using the same notation introduced before, this resulted in three data matrices 
each of dimension n = 618 and p = 26017. For each probe in each tissue, a linear 
regression model was fitted to regress out the effects of age and experimental batch, 
following the same procedure as in Grundberg et al. (2012). Residuals in adipose, 
LCL, and skin tissues were arranged in n x p matrices A^ 1 ), A^ 2 ), AT®, respectively, 
for further analysis using the proposed multiple-view matrix factorisation method. 

4.2 Experimental results 

Non-sparse MVMF was initially fitted for all combination of parameter pairs ( d , r) in 
a grid. For each model fit, we computed the percentage of variance explained in each 
tissue. These are shown in the 3D bar charts presented in Appendix, Section E, Figure 
9. The percentages of variance explained varied between 25.2% (d = r = 1, LCL) and 
87.3% (d = r = 160, skin). The following analyses are based on the d = r = 3 setting, 
which explains at least 40% of expression variance across tissues. Given that there 
are more than 26000 probes, and this is much larger than the sample size, this choice 
of parameters offers a good balance between dimensionality reduction and retaining a 
large portion of total variance. Although two other combinations of ( d,r ), i.e. (2,4) 
and (4,2), also explain a similar amount of total variance, we have found that the 
gene ranking results are not extremely sensitive to these values. For more details on 
this sensitivity analysis, see Appendix, Section C. 

The sparse version of our model, sMVMF, to each subsample in stability selection 
procedure to rank gene expressions explaining a large amount of shared and tissue- 
specific variances respectively. A detailed description of the procedure is presented 
in Section 2.3. In summary, 1000 random subsamples were generated each consisting 
of 309 subjects randomly and independently sampled without replacement from a 
total of 618. No twin pair was included in any subsample in order to remove pos¬ 
sible correlations due to zygosity. sMVMF was fitted to each subsample, where the 
sparsity parameters were fixed such that each column of the transformation matrices 
comprised exactly 100 non-zero entries. There were 3274 mRNA expression probes 
that were selected at least once from any of the transformation matrices. 

Probes that explain a large amount of expression variance exclusively in one tissue 
are of particular interest. To make such probes visually discernible we propose a new 
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SPOW Plot for TwinsUK Study 



Figure 4: TwinsUK study: resulting SPOW plot. The wheel comprises four rings, 
which correspond to shared, adipose-, LCL-, and skin-specific variability from the in¬ 
ner ring. It is also evenly divided into 3274 fan slices, corresponding to 3274 mRNA 
expression probes that are selected at least once in all subsamples. Probes are re¬ 
ordered by their selection probabilities in the transformation matrix in the shared 
component. Brighter colour denotes higher probability, whereas darker colour de¬ 
notes lower probability. We are particularly interested in probes with high selection 
probability exclusively in one ring. 


visualisation tool, the SPOW (Selection Probability Wheel) plot. The plot in Figure 
4 consists of 3274 fan slices corresponding to probes that are selected at least once 
in all subsamples, re-ordered by their selection probabilities in V*. The wheel is 
further divided into four rings, representing shared, adipose-, LCL-, and skin tissue, 
respectively. Each ring is assigned a unique colour spectrum to illustrate selection 
probabilities of the probes: brighter colours denote a higher probability and darker 
colours denote a lower probability. Probes featuring exclusively shared or tissue- 
specific variability can be found along the radii where only one part is painted in a 
bright colour and the other three parts are colored in black. The SPOW plots for 
the top 200 probes that explain shared and tissue-specific variability respectively are 
presented in Appendix, Section E, Figures 10 to 13, where such probes can be more 
easily captured. 
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Four groups of mRNA expressions were selected for further investigation, corre¬ 
sponding to shared-exclusive, adipose-, LCL-, and skin-exclusive expressions. Each 
group consisted of probes whose selection probabilities were larger than 0.5 in the 
corresponding transformation matrix and less than 0.005 in the other transformation 
matrices. These thresholds were set to give a manageable number of featured gene 
probes while tolerating occasional selection in the other groups. This procedure se¬ 
lected 294 genes for further study, including 114 adipose-exclusive, 83 LCL-exclusive, 
64 skin-exclusive, and 33 shared-exclusive genes. We summarise the results in Table 
1. A Venn-diagram representation of the results is given in Appendix, Section D. 

Table 1: TwinsUK study: summary of results. There are additionally 33 shared- 
exclusive genes. 



% of variance 
explained by 
tissue-specific 
component 

% of variance 
explained by 
shared 
component 

Number 
of tissue- 
exclusive 
probes 

Number 
of tissue- 
exclusive 

genes 

Adipose 

27.0 

14.7 

132 

114 

LCL 

30.8 

12.1 

91 

83 

Skin 

32.6 

11.5 

74 

64 


For each tissue, we performed an enrichment test by overlapping genes in our 
list with genes contained in the TiGER and VeryGene databases to examine the 
extent of agreement. In addition, a Gene Ontology (GO) biological process pathway 
enrichment test (Ashburner et al, 2000) and a Cytoscape pathway (CP) analysis 
(Saito et al., 2012) were carried out to reveal the function of the pathways which the 
261 tissue-exclusive genes belonged to, and FDR-corrected p-values were reported 
(See Supplementary Material, Table T1 and T2 for full results). Below we present 
test results for each group of genes separately for each tissue. We also report the 
selection probability (SP) for some selected probes. 

Skin-exclusive genes. 

15 of the 64 genes from our skin-exclusive list are contained in the combined TiGER/VeryGene 
list, giving rise to significant enrichment of our list with Fisher exact test p-value 
p < 10~ 16 . The overlapping genes include serine protease family genes KLK5 (SP: 

1.000) and KLK7 (SP: 1.000), which are highly expressed in the epidermis and re¬ 
lated to various skin conditions, such as cell shedding (desquamation) (Brattsand 
and Egelrud, 1999). Another member ALOX12B (SP: 1.000) controls producing 
12R-LOX, which adds an oxygen molecule to a fatty acid to produce the 12R- 
hydroperoxyeicosatetraenoic acid that has major function in the skin cell proliferation 
and differentiation (de Juanes et al, 2009). The skin-exclusive genes have also been 
found significantly enriched in two biological processes, namely epidermis develop¬ 
ment and cell-cell adhesion (p < 0.001 and p = 0.03, respectively). 
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LCL-exclusive genes. 

LCLs are not natural human cells: they are laboratory induced immortal cells that 
have abnormal telomerase activity and tumorigenic property (Sie et al., 2009). Since 
neither TiGER nor VeryGene assessed transcriptomic profile in LCL cells, we obtained 
LCLs data from Li et al. (2010), in which the authors compared LCLs expression 
profile in four human populations and reported 282 LCL specific expression genes. 
9 of those genes are contained in our LCL-exclusive gene list, giving a Fisher exact 
test p < 10 -16 . These include CDK5R1 (SP: 0.961) and HEY1 (SP: 1.000), which 
are key genes in the transformation of B lymphocytes to LCLs (Zhao et al., 2006). 
Pathway analysis of the LCL-exclusive genes reveals several aging and cell-death 
related pathways such as regulation of telomerase (CP enrichment test, p = 0.014), 
small cell lung cancer (CP enrichment test, p = 0.019), and cell cycle checkpoints 
(CP enrichment test, p = 0.021). These results show that our tissue-exclusive genes 
represent tissue unique molecular functions and biological pathways, which may be 
used to validate known pathways or discover new biological mechanisms. 

Adipose-exclusive genes. 

ApoB (SP: 1.000) is the only member in our adipose-exclusive list which is also 
contained in the list of known adipose-specific expression genes (Fisher exact test, 
p = 0.05). ApoB is one of the primary apolipoproteins that transport cholesterol to 
peripheral tissues (Knott et al., 1986) and it has been widely linked to fat formation 
(Riches et al., 1999). In adipose, the selected genes are found significantly enriched in 
triglyceride catabolic process pathway (p = 0.022), which is in line with the fact that 
adipose tissue is the major storage site for fat in the form of triglycerides. Pathway 
analysis reveals that genes in the adipose-exclusive list are significantly enriched in 
triglyceride catabolic process pathway (p = 0.022), which agrees with the fact that 
adipose tissue is the major storage site for fat in the form of triglycerides. In addition, 
these genes are enriched in inflammation pathways, such as lymphocyte chemotaxis 
(p = 0.016) and neutrophil chemotaxis (p = 0.027). This coincides with previous 
findings of the complex and strong link between metabolism and immune system in 
adipose tissue (Tilg and Moschen, 2006). 

For this tissue we were also able to further investigate the causes for the observed 
adipose-exclusive gene expression variability. One possible explanation could be that 
environmental factors influenced an individual’s epigenetic status, which subsequently 
regulated gene expression (Razin and Cedar, 1991). As a mediator of gene regulatory 
mechanisms, DNA methylation is crucial to genomic functions such as transcrip¬ 
tion, chromosomal stability, imprinting, and X-chromosome inactivation (Lokk et al., 
2014), which consequently influence an individual’s tissue development (Ziller et al., 
2013). It thus seemed reasonable to hypothesise that the expression of tissue-exclusive 
genes could be modified by their methylation status in the same tissue. 

We sought to identify genes featuring a statistically significant linear relationship 
between the gene’s methylation profile and its expression value from the same tissue. 
In adipose biopsies, where both transcriptome and methylation data is available, we 
found that 68.4% (78 out of 114 genes) of the genes had expression levels significantly 
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associated with their methylation status using a linear fit (Bonferroni correction, 
p < 0.05) (See Supplementary Material, Table T3, for full lists). We then wanted to 
assess whether a similar number of linear associations could be found by chance only 
by randomly selecting any genes, not only those that feature adipose-exclusive vari¬ 
ability, and testing for association between gene expression and methylation levels. 
This was done by randomly extracting the same, fixed number (132) of expression 
probes and corresponding methylation levels from adipose tissue, and fitting a linear 
model as before. By repeating this experiment 1000 times, we obtained the empirical 
distribution reported in Appendix, Section E, Figure 14. This distribution suggested 
that all the proportions were below 0.2, compared to our observed proportion of 
0.684, which provided overwhelming evidence that DNA methylation was an impor¬ 
tant factor affecting the expression of the tissue-exclusive genes. It was notable that 
the adipose-exclusive variability of ApoB was regulated by methylation at 50bp up¬ 
stream of the Transcriptional Starting Site (linear fit, p = 2.1 x 10~' 5 ), which agreed 
with the findings that the promoter of ApoB has tissue-specific and species-specific 
methylation property (Apostel et al, 2002). Apart from ApoB, we also found that 
methylation in Syk was associated with Syk expression level, which was potentially 
involved in B cell development and cell apoptosis (Ma et al., 2010). 


5 Conclusion and Discussion 

The proposed sMVMF method facilitates the comparison of gene expression vari¬ 
ances across multiple tissues. The primary challenge of this task arises from the 
interference between substantial co-variability of gene expressions across all tissues 
and substantial variability of gene expressions featured only in specific tissues. Char¬ 
acterising tissue-specific variability can shed light on the biological processes involved 
with tissue differentiation. Analysing shared variability can potentially reveal genes 
that are involved in complex or basic biological processes, and may as well enhance 
the estimation of tissue-specific variability. 

sMVMF has been used here to compare gene expression variances in three human 
tissues from the TwinsUK cohort. 261 genes having substantial expression variabil¬ 
ity exclusively featured in one tissue have been identified. Enrichment tests showed 
significant overlaps between our lists of tissue-exclusive genes and those reported in 
the TiGER and VeryGene databases, which were established by comparing mean ex¬ 
pression levels. This confirms the link between tissue-specific expression variance and 
the biological functions associated with particular tissues. In future work, it would 
be interesting to explore the functions of the tissue-exclusive genes from our list that 
have not been reported in existing databases. We further showed adipose-exclusive 
expression variability was driven by an epigenetic effect. Using these results as a 
guiding principle, we expect our methods and results could improve efficiencies in 
mapping functional genes by reducing the multiple testing and enhancing the knowl¬ 
edge of gene function in tissue development and disease phenotypes. Future works 
would consist of investigating the outcome of tissue-exclusive expression variability, 
for which we can perform association studies between expressions of tissue-exclusive 
genes and disease phenotypes related to adipose and skin tissues. 
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Appendix A Simulation setting 

As introduced in Section 3.1 of the main text, variance of the first 30 variables 
(columns) in the random matrices (m = 1,2,3) are controlled by three latent 

factors: Hi, H 2 , H 3 , which are real valued univariate random variables generated 
from independent normal distributions as follows: 

Ri~AA(0,5 2 ) ; H 2 ~ A/"(0,3.5 2 ) ; # 3 ~ AA( 0 , 2 2 ) (19) 

where cr 2 ) refers to normal distribution with mean /r and standard deviation a. 

Variance of the first 30 variables (columns) in the random matrices Z (m = 
1,2,3) is controlled by three latent factors: hi, h 2 , h 3 , where h m only affects Z^ m \ 
These latent factors are also generated from independent normal distributions: 

hi ~ A7(0,2.8 2 ) ; h 2 ~ Af(0,3.2 2 ) ; h 3 ~AA(0,3 2 ) (20) 

The latent variables in (19) and (20) control the variance of the first 30 variables in Y 
and Z via some constant factors which we shall define. Specifically, each value in the 
first 30 columns of Y is obtained by multiplying one latent variable from {H 1 , H 2 , i7 3 } 
with a constant factor from one of the two row vectors a or (3, so that the variance 
pattern in Y^ m > is precisely as is illustrated in the middle panel in Figure 2 of the main 
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paper. Similarly, each value in the first 30 columns of Z is obtained by multiplying 
one latent variable from {h\, h 2 , h 3 } with a constant factor from one of the row vectors 
7ij 72 ) 73 ) such that the variance pattern in Z is precisely as is illustrated in the 
right panel in Figure 2 of the main paper. The details are given as follows: 

a = (0.3, 0.5,0.6,0.8,1,1,0.8, 0.6,0.5,0.3) 

/3 = (0.6,0.7,0.8, 0.9,1,1, 0.9,0.8,0.7, 0.6) 

71 = (vi,vi,vi,vi,vi), where vi = (1,1/3, 2/3,1, 2/3,1/3) (21) 

72 = (v 2 ,v 2 ,v 2 ,v 2 ,v 2 ), where v 2 = (2/3,1,1/3,1/3,1, 2/3) 

13 = (v 3 ,v 3 ,v 3 ,v 3 ,v 3 ), where u 3 = (1/3, 2/3,1, 2/3,1/3,1) 


Let denote the ( i,j) th entry of Y^ m \ Our simulated data are generated as 

follows: for i = 1, 2,..., 100 and for m = 1,2, 3: 


1. Generate H\, H 2 , H 3 , hi, h 2 , and h 3 according to (19) and (20). 

2. Generate g 00 from independent normal distributions with zero mean and 
variance < 7 /, where < 7 / = 1 in setting I and a/ = 4 in setting II. 


3. Compute/Set: 


7”i’o = a ' Hi 
YlLo = a ■ H 2 

=/»•«> 

\d m ) _ rv 

1 i,31:500 — u 

y(m) _ 7 

■^1,1:30 — 3/m ' n-m 

y( m ) _ n 
^1,31:500 “ U 


Finally, compute: X( m ) = Y ("*) + Z^ . 


Appendix B Additional simulation 

In this additional simulation we use the same settings as in the previous section except 
that is generated from independent normal distributions with zero mean and 

variance 4. The same type of heatnraps as in Figure 2 and 3 of the main paper are 
produced and presented in Figure 5 and 6 respectively. We can visually conclude that 
sMVMF remains the best model in identifying the variables which drive shared and 
tissue-specific variance. Remarkably, Levene’s test hardly detects any genes whose 
variance is significantly larger than the corresponding genes in the other tissues due 
to increased noise level. 
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True and estimated latent structure controlling shared variance 


(A) shared LF structure 


LF l 
LF 2 
LF 3 



d = 1 
d = 2 
d = 3 


Active Inactive 


(B) sMVMF 



Top 10 
Top 20 
Top 30 


(C) Stacked-PCA 

P T 
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Figure 5: Each latent factor (LF) is only active in a block of 10 signal-carrying vari¬ 
ables, and controls the amount of variance of those variables that is shared amongst 
datasets. The (A) panel shows the true latent structure used to generate the data. 
Panels (B) and (C) show the estimated probabilities that each variable has been 
selected as signal-carrier using sMVMF and a stacked-PCA approach, respectively. 
sMVMF best captures the true shared LF structure whereas stacked-PCA tends to 
identify variables with large variance but fails to identify the LF structure. 


Selected variables controlling tissue-specific variance 

Tissue 1 
Tissue 2 
Tissue 3 


Figure 6: Three different methods - sMVMF, Levene’s test and PCA - are used 
to detect random variables whose variance pattern is dataset-specific. Each heatmap 
represents the selection probabilities estimated by each method: (A) sMVMF produces 
patterns that best match the true tissue-specific variances shown in the right panel of 
Figure 1 in the main paper; (B) Levene’s test hardly detects any gene whose variance is 
significantly larger than the corresponding genes in the other tissues due to increased 
noise level; (C) The PCA-based method cannot distinguish between shared and tissue- 
specific variability, and fails to recover the true pattern on variables with large shared 
variance. 



Appendix C Robustness study on the choice of 

(d, r) 

To investigate the robustness of (d, r ) on selected (shared- and tissue- exclusive) 
genes would require re-running the full analysis on all (d, r ) pairs on the 11x11 grid 
considered in our analysis, which would involve very intensive computation. Here we 
present a study in smaller scale in which we restrict the total amount of variance 
explained in adipose tissue to about 42%, and this gives us three pairs of (d, r ): (3, 3) 
which was the pair used to fit the sMVMF to identify shared- and tissue- exclusive 
genes in the paper, (2,4) and (4, 2). We present the percentages of shared and tissue- 
specific variance explained for these three combinations of (d,r) in Table 2. The 
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figures in LCL and skin tissues are very similar (±3%) to the adipose tissue for each 
(d,r) 


Table 2: Percentage of variance explained in adipose tissue 


(d,r) 

By shared component 

By tissue-specific component 

Total 

(3,3) 

14.7 

27.0 

41.7 

(2,4) 

10.3 

32.6 

42.9 

(4,2) 

23.3 

18.4 

41.7 


Notably, although the percentages of explained variance are approximately equal 
for the three combinations considered, the percentages within each component (shared 
and tissue-specific variances) vary substantially, in particular between (2,4) and (4, 2). 
Therefore, this incomplete comparison seems to give a valid illustration of the robust¬ 
ness of gene selection results on the full grid of ( d,r ). 

To evaluate the robustness of the shared- and tissue- exclusive genes with respect 
to the choice of ( d,r ), we repeated our analysis for (d, r) = (2,4) and (d, r) = (4,2) 
in the same way as for (d, r ) = (3, 3). We adjusted the selection criteria (threshold of 
selection probabilities) following the same principle as introduced in the paper so that 
the same number of shared- and tissue- exclusive genes (±1 when there were ties) were 
selected as in the lists for (d, r) = (3, 3). We present the Venn diagrams summarising 
the overlaps between the three combinations of (d, r) parameters in Figure 7. 

The results showed that adipose- and skin- exclusive genes were very robust to 
the choice of (d, r) since more than 77% of genes appeared in the lists obtained from 
all three combinations of (d, r). The shared-exclusive genes were fairly robust to 
the choice of (d, r) in that there were more than 90% of overlaps between the lists 
obtained from (d, r) = (3,3) and (d, r) = (2,4), and about 40% of overlaps with 
the list obtained from (d, r) = (4,2). However, the percentage of overlaps would 
increase to 70% if we restrain the comparison to the top 20 shared-exclusive genes. 
For LCL-exclusive genes, there were about 40% of overlaps among the three pairs of 
(d, r). Moreover, the highlighted genes mentioned in the main paper were all retained 
in the lists of genes selected using the other combinations of (d, r), except for the 
LCL-exclusive gene CDK5R1 which was absent from (d,r) = (4,2). We therefore 
conclude that given the substantial difference in the percentages of shared and tissue- 
specific variance explained using different combinations of (d, r), the lists of shared- 
and tissue- exclusive genes were robust to the choice of (d, r), in particular if such 
lists were small. 


Appendix D Venn-diagram analysis 

We present a Venn diagram in Figure 8 summarising our findings from the TwinsUK 
analysis. As mentioned in the main paper, we identified 114 adipose-exclusive, 83 
LCL-exclusive, and 64 skin-exclusive genes. In addition, 33 genes which drove the 
shared variability across all three tissues yet without driving tissue-specific variability 
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Shared-excl. genes 


adipose-excl. genes 



LCL-excl. genes 





d = r = 3 d = 2, r = 4 d = 4, r = 2 


Figure 7: The Venn diagrams summarise the overlaps of the shared- and tissue- 
exclusive genes selected for three combinations of (d, r ) pairs which explain about 
42% of the total variance in the adipose tissue. These plots showed that adipose- and 
skin- exclusive genes were very robust to the choice of (d, r) and shared- and LCL- 
exclusive genes were fairly robust. 


in any tissue were identified. Moreover, 2 genes (“AQP9” and “TYMP”) were identified 
to have driven adipose- and LCL-specific variance but not skin-specific variance, while 
4 genes (“CCND1”, “GPC4”, “GSDMB”, and “TUBB2B”) were found to have driven 
adipose- and skin-specific variance but not LCL-specific variance. 


26 

















Venn-diagram of TwinsUK analysis 



Figure 8: The Venn diagram shows that there were 114 adipose-exclusive, 83 LCL- 
exclusive, 64 skin-exclusive, and 33 shared-exclusive genes extracted from our analysis. 
Using the SPOW plot in Figure 4 of the main paper, we were also able to identify 
2 genes (“AQP9” and “TYMP”) which drove tissue-specific variability in adipose and 
LCL tissues but not in skin; in addition we also identified 4 genes (“CCND1”, “GPC4”, 
“GSDMB”, and “TUBB2B”) which drove tissue-specific variability in adipose and skin 
tissues but not in LCL. 


Appendix E Plots 
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Figure 9: TwinsUK study: 3D boxplot showing the percentage of expression variance 
explained in adipose, LCL, and skin tissues on a grid of (d, r ) using the non-sparse 
MVMF. d is the total number of PPJs in the shared variance component, and r is the 
total number of PPJs in the tissue-specific variance component. The percentages vary 
between 25.2% (d = r = 1, LCL) and 87.3% (d = r = 160, skin). 
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shared adipose LCL skin 

Figure 10: TwinsUK study: SPOW plot (d = r = 3). The wheel contains the top 
200 most frequently selected probes from the transformation matrix in the shared 
component. We extract probes with bright colour in the shared variability (green) 
ring and dark colours in the other rings. 


28 



SPOW Plot - Top 200 probes in adipose-specific list 



shared adipose LCL skin 

Figure 11: TwinsUK study: SPOW plot (d = r = 3). The wheel contains the top 200 
most frequently selected probes from the transformation matrix in the adipose-specific 
component using sMVMF. We extract probes with bright colour in the adipose-specific 
(yellow) ring and dark colours in the other rings. 
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SPOW Plot - Top 200 probes in LCL-specific list 



shared adipose LCL skin 

Figure 12: TwinsUK study: SPOW plot (d = r = 3). The wheel contains the top 200 
most frequently selected probes from the transformation matrix in the LCL-specific 
component using sMVMF. We extract probes with bright colour in the LCL-specific 
(purple) ring and dark colours in the other rings. 
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SPOW Plot - Top 200 probes in skin-specific list 



shared adipose LCL skin 

Figure 13: TwinsUK study: SPOW plot (d = r = 3). The wheel contains the top 
200 most frequently selected probes from the transformation matrix in the skin-specific 
component using sMVMF. We extract probes with bright colour in the skin-specific 
(cyan) ring and dark colours in the other rings. 
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Figure 14: Proportion of randomly chosen genes for which the corresponding gene 
expression shows a significant linear association with the methylation probe. The 
experiment consists of 1000 random draws, and each draw involves 132 randomly cho¬ 
sen expression probes, which are tested for linear association with the corresponding 
methylation profiles. We conclude that observing a proportion as large or larger than 
0.684, which is what we obtained for our adipose-exclusive genes, is unlikely to happen 
by chance only. 
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